home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / srch.c < prev    next >
C/C++ Source or Header  |  1992-05-24  |  8KB  |  352 lines
  1. /* this file contains functions to support iterative ephem searches.
  2.  * we support several kinds of searching and solving algorithms.
  3.  * values used in the evaluations come from the field logging flog.c system.
  4.  * the expressions being evaluated are compiled and executed from compiler.c.
  5.  */
  6.  
  7. #include <stdio.h>
  8. #include <math.h>
  9. #include "screen.h"
  10.  
  11. extern char *strcpy();
  12.  
  13. static int (*srch_f)();
  14. static int srch_tmscalled;
  15. static char expbuf[NC];        /* [0] == '\0' when expression is invalid */
  16. static double tmlimit = 1./60.;    /* search accuracy, in hrs; def is one minute */
  17.  
  18.  
  19. srch_setup()
  20. {
  21.     int srch_minmax(), srch_solve0(), srch_binary();
  22.     static char *chcs[] = {
  23.         "Find extreme", "Find 0", "Binary", "New function", "Accuracy",
  24.         "Stop"
  25.     };
  26.     static int fn;    /* start with 0, then remember for next time */
  27.  
  28.     /* let op select algorithm, edit, set accuracy
  29.      * or stop if currently searching
  30.      * algorithms require a function.
  31.      */
  32.     ask:
  33.     switch (popup(chcs, fn, srch_f ? 6 : 5)) {
  34.     case 0: fn = 0;
  35.         if (expbuf[0] == '\0')
  36.             set_function();
  37.         srch_f = expbuf[0] ? srch_minmax : (int (*)())0;
  38.         if (srch_f)
  39.             break;
  40.         else
  41.             goto ask;
  42.     case 1: fn = 1;
  43.         if (expbuf[0] == '\0')
  44.             set_function();
  45.         srch_f = expbuf[0] ? srch_solve0 : (int (*)())0;
  46.         if (srch_f)
  47.             break;
  48.         else
  49.             goto ask;
  50.     case 2: fn = 2;
  51.         if (expbuf[0] == '\0')
  52.             set_function();
  53.         srch_f = expbuf[0] ? srch_binary : (int (*)())0;
  54.         if (srch_f)
  55.             break;
  56.         else
  57.             goto ask;
  58.     case 3: fn = 3; srch_f = 0; set_function(); goto ask;
  59.     case 4: fn = 4; srch_f = 0; set_accuracy(); goto ask;
  60.     case 5: srch_f = 0; srch_prstate(0); return;
  61.     default: return;
  62.     }
  63.  
  64.     /* new search */
  65.     srch_tmscalled = 0;
  66.     srch_prstate (0);
  67. }
  68.  
  69. /* if searching is in effect call the search type function.
  70.  * it might modify *tmincp according to where it next wants to eval.
  71.  * (remember tminc is in hours, not days).
  72.  * if searching ends for any reason it is also turned off.
  73.  * also, flog the new value.
  74.  * return 0 if caller can continue or -1 if it is time to stop.
  75.  */
  76. srch_eval(mjd, tmincp)
  77. double mjd;
  78. double *tmincp;
  79. {
  80.     char errbuf[128];
  81.     int s;
  82.     double v;
  83.  
  84.     if (!srch_f)
  85.         return (0);
  86.  
  87.     if (execute_expr (&v, errbuf) < 0) {
  88.         srch_f = 0;
  89.         f_msg (errbuf);
  90.     } else {
  91.         s = (*srch_f)(mjd, v, tmincp);
  92.         if (s < 0)
  93.         srch_f = 0;
  94.         (void) flog_log (R_SRCH, C_SRCH, v, "");
  95.         srch_tmscalled++;
  96.     }
  97.  
  98.     srch_prstate (0);
  99.     return (s);
  100. }
  101.  
  102. /* print state of searching. */
  103. srch_prstate (force)
  104. int force;
  105. {
  106.     int srch_minmax(), srch_solve0(), srch_binary();
  107.     static (*last)();
  108.  
  109.     if (force || srch_f != last) {
  110.         f_string (R_SRCH, C_SRCHV,
  111.             srch_f == srch_minmax   ? "Extrema" :
  112.             srch_f == srch_solve0   ? " Find 0" :
  113.             srch_f == srch_binary ?   " Binary" :
  114.                           "    off");
  115.         last = srch_f;
  116.     }
  117. }
  118.  
  119. srch_ison()
  120. {
  121.     return (srch_f != 0);
  122. }
  123.  
  124. /* display current expression. then if type in at least one char make it the
  125.  * current expression IF it compiles ok.
  126.  * TODO: editing?
  127.  */
  128. static
  129. set_function()
  130. {
  131.     static char prompt[] = "Function: ";
  132.     char newexp[NC];
  133.     int s;
  134.  
  135.     f_prompt (prompt);
  136. #ifdef AMIGA
  137.     cwrite(expbuf);
  138. #else
  139.     (void) fputs (expbuf, stdout);
  140. #endif
  141.     c_pos (R_PROMPT, sizeof(prompt));
  142.  
  143.     s = read_line (newexp, PW-sizeof(prompt));
  144.     if (s >= 0) {
  145.         char errbuf[NC];
  146.         if (s > 0 && compile_expr (newexp, errbuf) < 0)
  147.         f_msg (errbuf);
  148.         else
  149.         (void) strcpy (expbuf, newexp);
  150.     }
  151. }
  152.  
  153. static
  154. set_accuracy()
  155. {
  156.     static char p[] = "Desired accuracy (         hrs): ";
  157.     int hrs, mins, secs;
  158.     char buf[NC];
  159.  
  160.     f_prompt (p);
  161.     f_time (R_PROMPT, C_PROMPT+18, tmlimit); /* place in blank spot */
  162.     c_pos (R_PROMPT, sizeof(p));
  163.     if (read_line (buf, PW-sizeof(p)) > 0) {
  164.         f_dec_sexsign (tmlimit, &hrs, &mins, &secs);
  165.         f_sscansex (buf, &hrs, &mins, &secs);
  166.         sex_dec (hrs, mins, secs, &tmlimit);
  167.     }
  168. }
  169.  
  170. /* use successive paraboloidal fits to find when expression is at a
  171.  * local minimum or maximum.
  172.  */
  173. static
  174. srch_minmax(mjd, v, tmincp)
  175. double mjd;
  176. double v;
  177. double *tmincp;
  178. {
  179.     static double base;        /* for better stability */
  180.     static double x_1, x_2, x_3;    /* keep in increasing order */
  181.     static double y_1, y_2, y_3;
  182.     double xm, a, b;
  183.  
  184.     if (srch_tmscalled == 0) {
  185.         base = mjd;
  186.         x_1 = 0.0;
  187.         y_1 = v;
  188.         return (0);
  189.     }
  190.     mjd -= base;
  191.     if (srch_tmscalled == 1) {
  192.         /* put in one of first two slots */
  193.         if (mjd < x_1) {
  194.             x_2 = x_1;  y_2 = y_1;
  195.         x_1 = mjd; y_1 = v;
  196.         } else {
  197.         x_2 = mjd; y_2 = v;
  198.         }
  199.         return (0);
  200.     }
  201.     if (srch_tmscalled == 2 || fabs(mjd - x_1) < fabs(mjd - x_3)) {
  202.         /* closer to x_1 so discard x_3.
  203.          * or if it's our third value we know to "discard" x_3.
  204.          */
  205.         if (mjd > x_2) {
  206.         x_3 = mjd; y_3 = v;
  207.         } else {
  208.         x_3 = x_2;  y_3 = y_2;
  209.         if (mjd > x_1) {
  210.             x_2 = mjd; y_2 = v;
  211.         } else {
  212.             x_2 = x_1;  y_2 = y_1;
  213.             x_1 = mjd; y_1 = v;
  214.         }
  215.         }
  216.         if (srch_tmscalled == 2)
  217.         return (0);
  218.     } else {
  219.         /* closer to x_3 so discard x_1 */
  220.         if (mjd < x_2) {
  221.         x_1 = mjd;  y_1 = v;
  222.         } else {
  223.         x_1 =  x_2;  y_1 = y_2;
  224.         if (mjd < x_3) {
  225.             x_2 = mjd; y_2 = v;
  226.         } else {
  227.             x_2 =  x_3; y_2 = y_3;
  228.             x_3 = mjd; y_3 = v;
  229.         }
  230.         }
  231.     }
  232.  
  233. #ifdef TRACEMM
  234.     { char buf[NC];
  235.       sprintf (buf, "x_1=%g y_1=%g x_2=%g y_2=%g x_3=%g y_3=%g",
  236.                         x_1, y_1, x_2, y_2, x_3, y_3);
  237.       f_msg (buf);
  238.     }
  239. #endif
  240.     a = y_1*(x_2-x_3) - y_2*(x_1-x_3) + y_3*(x_1-x_2);
  241.     if (fabs(a) < 1e-10) {
  242.         /* near-0 zero denominator, ie, curve is pretty flat here,
  243.          * so assume we are done enough.
  244.          * signal this by forcing a 0 tminc.
  245.          */
  246.         *tmincp = 0.0;
  247.         return (-1);
  248.     }
  249.     b = (x_1*x_1)*(y_2-y_3) - (x_2*x_2)*(y_1-y_3) + (x_3*x_3)*(y_1-y_2);
  250.     xm = -b/(2.0*a);
  251.     *tmincp = (xm - mjd)*24.0;
  252.     return (fabs (*tmincp) < tmlimit ? -1 : 0);
  253. }
  254.  
  255. /* use secant method to solve for time when expression passes through 0.
  256.  */
  257. static
  258. srch_solve0(mjd, v, tmincp)
  259. double mjd;
  260. double v;
  261. double *tmincp;
  262. {
  263.     static double x0, x_1;    /* x(n-1) and x(n) */
  264.     static double y_0, y_1;    /* y(n-1) and y(n) */
  265.     double x_2;        /* x(n+1) */
  266.     double df;        /* y(n) - y(n-1) */
  267.  
  268.     switch (srch_tmscalled) {
  269.     case 0: x0 = mjd; y_0 = v; return(0);
  270.     case 1: x_1 = mjd; y_1 = v; break;
  271.     default: x0 = x_1; y_0 = y_1; x_1 = mjd; y_1 = v; break;
  272.     }
  273.  
  274.     df = y_1 - y_0;
  275.     if (fabs(df) < 1e-10) {
  276.         /* near-0 zero denominator, ie, curve is pretty flat here,
  277.          * so assume we are done enough.
  278.          * signal this by forcing a 0 tminc.
  279.          */
  280.         *tmincp = 0.0;
  281.         return (-1);
  282.     }
  283.     x_2 = x_1 - y_1*(x_1-x0)/df;
  284.     *tmincp = (x_2 - mjd)*24.0;
  285.     return (fabs (*tmincp) < tmlimit ? -1 : 0);
  286. }
  287.  
  288. /* binary search for time when expression changes from its initial state.
  289.  * if the change is outside the initial tminc range, then keep searching in that
  290.  *    direction by tminc first before starting to divide down.
  291.  */
  292. static
  293. srch_binary(mjd, v, tmincp)
  294. double mjd;
  295. double v;
  296. double *tmincp;
  297. {
  298.     static double lb, ub;        /* lower and upper bound */
  299.     static int initial_state;
  300.     int this_state = v >= 0.5;
  301.  
  302. #define    FLUNDEF    -9e10
  303.  
  304.     if (srch_tmscalled == 0) {
  305.         if (*tmincp >= 0.0) {
  306.         /* going forwards in time so first mjd is lb and no ub yet */
  307.         lb = mjd;
  308.         ub = FLUNDEF;
  309.         } else {
  310.         /* going backwards in time so first mjd is ub and no lb yet */
  311.         ub = mjd;
  312.         lb = FLUNDEF;
  313.         }
  314.         initial_state = this_state;
  315.         return (0);
  316.     }
  317.  
  318.     if (ub != FLUNDEF && lb != FLUNDEF) {
  319.         if (this_state == initial_state)
  320.         lb = mjd;
  321.         else
  322.         ub = mjd;
  323.         *tmincp = ((lb + ub)/2.0 - mjd)*24.0;
  324. #ifdef TRACEBIN
  325.         { char buf[NC];
  326.           sprintf (buf, "lb=%g ub=%g tminc=%g mjd=%g is=%d ts=%d",
  327.                 lb, ub, *tmincp, mjd, initial_state, this_state);
  328.           f_msg (buf);
  329.         }
  330. #endif
  331.         /* signal to stop if asking for time change less than TMLIMIT */
  332.         return (fabs (*tmincp) < tmlimit ? -1 : 0);
  333.     } else if (this_state != initial_state) {
  334.         /* gone past; turn around half way */
  335.         if (*tmincp >= 0.0)
  336.         ub = mjd;
  337.         else
  338.         lb = mjd;
  339.         *tmincp /= -2.0;
  340.         return (0);
  341.     } else {
  342.         /* just keep going, looking for first state change but we keep
  343.          * learning the lower (or upper, if going backwards) bound.
  344.          */
  345.         if (*tmincp >= 0.0)
  346.         lb = mjd;
  347.         else
  348.         ub = mjd;
  349.         return (0);
  350.     }
  351. }
  352.